A Renormalisation-Group Algorithm for Eigenvalue Density Functions of Interacting Quantum 

Systems 
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We present a certifiable algorithm to calculate the eigenvalue density function — the number of eigen- 
values within an infinitesimal interval — for an arbitrary ID interacting quantum spin system. Our method 
provides an arbitrarily accurate numerical representation for the smeared eigenvalue density function, 
which is the convolution of the eigenvalue density function with a gaussian of prespecified width. In 
addition, with our algorithm it is possible to investigate the density of states near the ground state. This 
can be used to numerically determine the size of the ground-state energy gap for the system to within a 
prespecified confidence interval. Our method exploits a finitely correlated state/matrix product state rep- 
resentation of the propagator and applies equally to disordered and critical interacting ID quantum spin 
systems. We illustrate our method by calculating an approximation to the eigenvalue density function for 
a random antiferromagnetic Heisenberg model. 

PACS numbers: 03.65.Bz, 89.70,+c 



The statics and dynamics of interacting quantum many- 
particle systems are still relatively poorly understood. In- 
deed, even calculating an approximation to such basic 
quantities as the ground-state energy appears to be ex- 
tremely difficult for many interesting systems. At least one 
reason for this is that for arbitrary local quantum systems 
this problem is complete for the complexity class QMA, 
which is the quantum analogue of NP ulSlj]- Of course, 
we do not really expect that there exist general efficient 
computational schemes to study eigenvalues and related 
thermodynamic properties. But, it is plausible that for re- 
alistic quantum systems there may exist efficient schemes 
to calculate certain physical properties like approximations 
to energy gaps and other thermodynamic properties. 

The development of the density matrix renormalisation 
group (DMRG) has provided us with what promises to 
be an efficient way to calculate physical properties of the 
ground states of interacting quantum systems in ID (see 
0| and references therein for a detailed description of the 
DMRG). While the method was originally developed to 
obtain approximations to the ground state of a regular in- 
teracting quantum spin lattice system in ID, the DMRG 
is an extremely flexible method and has been recently ex- 
tended to apply to a diverse number of situations, such as 
the calculation of short-time dynamics |0, ^.dissipation 
@h eigenstates with definite momentum |2|, and, re- 
cently, higher dimensions | ljj]. 

Whether the DMRG and related algorithms actually 
compute approximations to the ground state of a quan- 
tum system instead of low-lying excited states is an open 
question. This problem is difficult to answer because the 
DMRG cannot be certified, i.e., once the DMRG produces 
a ground-state approximation there is no way to prove that 
this approximation is correct to within some prespecificed 
confidence interval. However, this situation is changing; 
there have recently been several works which provide certi- 
fiable DMRG-like algorithms to approximate ground states 



of interacting spin systems fill 111 111 fLU 111 {iMvh. 

One situation where DMRG-related algorithms have 
been less successful is in the calculation of the eigenvalue 
density function /jLh{x) and the eigenvalue counting func- 
tion N H (x), which counts the number of eigenvalues of 
a hamiltonian H with value less than x 12111 . (The two 
functions are connected by Hh{ x ) = dN H {x)/dx.) Per- 
haps the closest general method which has been developed 
along these lines is due to Porras, Verstraete, and Cirac tJ]. 
This method calculates eigenstates with definite linear mo- 
mentum for ID quantum spin systems on a ring. However, 
it is possible for the method of Porras, Verstrate, and Cirac 
to miss eigenstates in the same way that the DMRG can 
miss the ground state and end up in a local minima. While 
this never appears to occur in practice it would be desir- 
able to have a method which trades this uncertainty against 
some other approximation. Additionally, if the method of 
Porras, Verstraete, and Cirac is used to calculate Hh(x) for 
x S> Oil) above the ground state energy then this would 
require exponential resources. 

In this Letter we introduce a certifiable method to cal- 
culate systematic approximations to /Xj? (a;) for rather gen- 
eral ID quantum spin systems (a generalisation to two and 
higher dimensions is available, which uses a slight modi- 
fication of the technology of Verstraete and Cirac fioll and 
fish ). We assume neither translation invariance nor non- 
criticality. Our method is an approximation because it cal- 
culates a "smeared" version Jl H {x) of Hn{x), which is the 
convolution of /ijj(ic) with a gaussian which has a width 
which can be reduced with a complexity that provably 
scales polynomially with n. Our method doesn't suffer 
from the uncertainty of the method of Porras, Verstraete, 
and Cirac, i.e., that maybe some eigenstates are missed. 
However, the price we pay for this is that while it is cer- 
tain that every eigenvalue is represented in fl H {x) there is 
some inevitable uncertainty in the calculated positions of 
the eigenvalues. 
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FIG. 1: The structure of the propagator e ltH for a quantum spin 
system for times t which are short in comparison to the size of 
the system (see II lL 1 1 9L 12011 for additional discussion). 

The outline of this Letter is as follows. We begin by in- 
troducing some definitions and we introduce the class of 
systems we study. We then describe our numerical renor- 
malisation group method to calculate Jl H (x). We conclude 
with some numerical results of our method applied to a 
random antiferromagnetic Heisenberg model. 

We will, for the sake of clarity, introduce and de- 
scribe our method for a chain of n distinguishable spin- 
| particles. Thus, the Hilbert space H for our system 

is given by H = (g)™ C 2 . Consider the C*-algebra 
B(TC) which is the Hilbert space of all (bounded) lin- 
ear operators A on 7i with inner product (A, B) = 
ti(A^B). An orfhonormal basis for B(7i) is given by 
by a a = a a ° <g) a ai •••cr Q —\ a 3 G Z/4Z, < 
j < n — 1, which we call the standard operator basis, 
wherea" = [(J ?) , (? J) , ($?) , (J _? x )] , is the vec- 
tor of Pauli sigma matrices. We write the structure con- 
stants g a/3 7 for the C*-algebra generated by a a : o a o~ 13 = 

X^7=o 5 ,Q/3 7 £j7 - The family H of local hamiltonians we 

focus on is defined by H = X^=o wnere hj 1S an 
interaction term which couples only neighbouring spins j 
and j + 1. We assume the standard energy normalisation 
whereby \\hj\\ = 0(1). The interaction hj may vary with 
position. 

The objective of this Letter is to understand the distri- 
bution of eigenvalues for the operator H. To do this we'll 
study the eigenvalue density function ^h(x) for H which 
is given by 

2»-l 

where 5(x) is the Dirac delta function and Ej are the eigen- 
values of H. The eigenvalue counting function Nh(x) of 
H is defined to be equal to the number of eigenvalues of H 
which are less than or equal to x. 

The eigenvalue density function for an opera- 

tor H has a delta function spike at the position of each 
eigenvalue of H. Notice that we have normalised the 
eigenvalue density function to have area 1, i.e. 

J_ dx ijlh{x) = 1. We have done this principally so that 
we can compare the eigenvalue densities for operators on 
different Hilbert spaces. The eigenvalue counting function 



N H (x) can be expressed in terms of the eigenvalue density 
function Hh(x) as N H (x) = 2" J_ duj [i h (uj). 

We obtain the eigenvalue density function (x) f° r an 
operator H via the following procedure. Write H in its 
eigenbasis, H = Ylj^ 1 Ej\Ef)(Ej\, and consider the 

propagator U (f) = e ltH = Y^Jq 1 ^ Eit \ E j){ E j\- Tak- 
ing the fourier transform of the scaled propagator ^U(t) 
yields 

±U(p) = ^[U{t)\ = ^Y / 6(E j -u)\E j )(E j \, 

3=0 

(2) 

where the fourier transform pair (#[•], 5 is defined 
to be F(u) = $[f(t)\ = f^dtfWe-™ and f(t) = 
$- x [F(uj)] = i f^dtF^e^. If we take the 
trace of ^U(uj) we find 2tt^l h (uj) = ti(U(uj)) = 

The calculations in the previous paragraph show that if 

we know ir(U(t)) for arbitrary times then we have enough 

information to extract /j>h(x). Now we show that if we 

only know an approximation V(t) to U(t) valid for some 

time \t\ < T then we can still extract an approximation 

Vn{x) to Hh(x) which can be systematically improved as 

T is increased. The idea is to introduce a scalar-valued 

windowing function Xt if) which cuts off the propagators 

U(t) and V{t) outside \t\ < T so that XT(t)V(t) ~ 

XT(t)U(t) for all t. One convenient choice for Xr(t), 

which we use in the sequel, is the gaussian: 

_ t 2 
e it 2 

If we now study ^Xr(t) tr(?7(t)), rather than 
t^- tr (C/(t)), then a fourier transform and an application 
of the convolution theorem yields 

1 27T 2 "- 1 

-fc(t)tr([/(t))] = _£ Xt(E 3 - u), (4) 

3=0 

where Xt(^) is the fourier transform of the window- 
ing/characteristic function. It is straightforward to identify 
Eq. © as a convolution 

^ff[XT(t)tr(E/(t))] = 2ir(xr*MH)(w). (5) 

In this way we identify the fourier transform of 
^tXt(^) tr(C/(t)) with a smearing of the eigenvalue den- 
sity function /j.h(u) with a smearing function Xt{^)- 
When Xr(t) is chosen to be a gaussian, as in Eq. 0, it is 
clear that increasing the time window T reduces the width 
of x t (uj). Thus, in the limit T — > oo we smoothly (but 
not uniformly!) recover the eigenvalue distribution func- 
tion: fi H (uj) = lim T ^ 00 (xT *//ir)(w). 

It is immediate that if we now have only an approxima- 
tion V(t) to U(t) which is good |H for \t\ < T then 
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the fourier transform of g(t) = ^Xr{t) tr(V(t)) will be 
close to that of f(t) = ^Xt(0 tr(Z7(i)). One way to jus- 
tify this is either to exploit the Parseval's relation, i.e., that 
the fourier transform is a unitary operation on L 2 , or to use 
the result that if ||/ — <?||i < e then || jf — fflloo — e > where 
|| • ||i and || • Hoc denote the L 1 and norms, respectively. 

We can also study the trace of the propagator in imagi- 
nary time: consider U(t + if3) = e - ^ e . Taking the 
trace yields 

2 n — 1 

/(it + 0) = tx(U(t + = e-W'e™'*. (6) 

We obtain /j, h from /(it + (3) for a fixed /3 by computing 
the laplace transform inversion integral 

/•/3+ioo 

Hb(") = / dsf(s)e s ". (7) 

•/ /3 — too 

For a fixed /3 this can be done with an inverse fourier trans- 
form in t: fi H (u}) = e /3w 5 r ~ 1 [/(it + /?)]. It is straight- 
forward to see that if we only know an approximation 
V(t + i(3)toU(t + i(3) valid for \t\ < T then the fourier 
transform g(u) along t of the cutoff trace g(it + /?) = 
XT(^)t r (^(i + provides a good approximation to 
e -0u(XT * Hh)[(jS). Thus, after normalisation, e^ u "g(uj) 
provides a good representation for the smeared eigenvalue 
density function for uj ~ E a , where E is the ground state 
energy. This representation becomes exponentially worse 
as u increases. By combining this representation with the 
one obtained from pure time evolution allows us to tradeoff 
the errors in the two representations. 

The preceding discussion serves to establish the fact that 
a good approximation V(t) to U(t) = e lHt which is valid 
for complex times \t\ < T provides sufficient information 
to resolve the eigenvalue distribution function on a length- 
scale 5 ~ 0(j:). In the next part of this Letter we pro- 
vide a numerical algorithm, closely related to the DMRG, 
which efficiently calculates a numerical representation for 

(Xt*Mh)H- 

The crucial idea underlying our numerical method is that 
a good approximation V(t) to the propagator U(t) for a 
local ID quantum spin lattice system can be stored effi- 
ciently (i.e. with polynomial resources in n) on a classical 
computer for for \t\ < T, where T ~ 0(log(n)) ifTTIl 
(see also 11191 |2(J]). See Fig. ^ f° r an illustration of the 
structure of the propagator e ltH for an arbitrary local spin 
system. This result allows us to certify that our algorithm 
can correctly obtain an approximate representation for the 
smeared eigenvalue distribution function to within a con- 
stant lengthscale which can be arbitrarily large (but scaling 
at most logarithmically with n.) 

The way we actually store a representation for V(t) is 
as a matrix product operator (MPO) ITila]. What we mean 
by this is that we represent an operator W G B(7i) in the 
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FIG. 2: The (scaled) eigenvalue density function for an in- 
stance of the random antiferromagnetic Heisenberg model H = 
J2k=o Jk&k ■ °"fc+i> where the couplings Jj, were chosen uni- 
formly at random from the interval (0, 1), for a chain of n — 20 
spins. The eigenvalue density function has been scaled so that the 
total area beneath the curve is 2 20 , the total number of eigenval- 
ues. The calculations were performed with maximum auxiliary 
dimension D — 20 and T = 20. The windowing gaussian has 
width s = 5, thus the eigenvalue density function is correct down 
to a lengthscale 27r/5. The theoretical worst-case error arising 
from the Trotter decomposition and truncation of D is O(10 -2 ) 
according to the L 2 and L>oo norms. 

following fashion 

W = V A j "A 31 • • • A^-V ®^' 1 ®- • -®o- jn -\ (8) 
jeQ„ 

where Q n = (Z/4Z) x,i , and A jo and A^'- 1 are a col- 
lection of four C x 1 sized row vectors (respectively, four 
1 x D n _i sized column vectors) and A Jfc are four x D k 
sized matrices. Note that Ck+i = D k . The dimensions 
Cfc and Dk are called the auxiliary dimensions for site 
k. It is clear that if the sizes of the auxiliary dimensions 
are bounded by polynomials in n, i.e. Ck < poly(n) and 
Dk < poly(n), then the operator W can be stored with 
polynomial resources in n. Also note that all operators can 
be represented exactly as in Eq. (|8j by taking the auxiliary 
dimensions to be large enough: C k = D k = 2™ suffices. 

We obtain the MPO representation for V(t) via the fol- 
lowing method. First we break H into two pieces A and 
B which contain the interaction terms on the even (respec- 
tively, odd) sites. Note that each term within A (respec- 
tively B) commutes with all the other terms within A (re- 
spectively, B). Then we exploit the Lie-Trotter expansion 

e lHt = lim (e lA ™e lB ™) m , (9) 

m — >oo 

to write our expression for V(t), i.e., we pick some m and 
write V(t) = W m , where W = e lA ™e iB ™. We next 
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FIG. 3: The Schmidt coefficients as a function of time between 
the first 10 spins and the second 10 spins for the MPO rep- 
resentation of the propagator e ltH for the random Heisenberg 
model H studied in Fig. |2| (The decay of the Schmidt co- 
efficients provide a measure of how difficult the propagator is 
to store 1 16] and are calculated as the squares of the singular 
values of M jk = A. jo A h ■ ■ ■ A j9 A k ° ■ ■ ■ A fcfl , where V(t) = 
Ej G Q 20 A j0 A J1 ■ ■ ■ A jl9 a j ° g> a n ® ■ ■ ■ <8 a 119 .) The Schmidt 
coefficients grow approximately exponentially as a function of 
time. 

express e lhl ™ = Y? a ,p=o c a/3°? ® a j+i- B ^ writing W in 
terms of the standard product operator basis we obtain the 
MPO representation 

W = BJ ° Ah ■ ■ ■ B^-V <g> a jl (g> ■ ■ ■ ® a 3n - x , 
jeQ n 

(10) 

where Bi° = c[- „, A Jk = s jk c^ on even sites, A Jk = 
gj* c [fe] on odd s jt eS; an( j B^- 1 = 8j n _ ua . 

Now we show that if two MPO's J and K have 
maximum auxiliary dimensions Dj and .D#- then 
Jif is expressible as a MPO L with auxiliary di- 
mension D JK < DjD K . Representing J = 
J2jeq n Ai ° AJl • • • A 3 "- 1 a 30 (g) cr J1 g> • • • <g> o^"- 1 and 
K - Ej G Q„ B^B^' 1 • • • B jn - 1 a jo (g> a jl (g> ■ ■ ■ <8> 
Q-Jn-i an( j taking the product gives L = JK = 
SjeQ„ C jo C jl ■ ■ ■ C jn - l a jo (g> o n <g> ■ ■ ■ <g> a 3 — 1 , where 

3 

c l -< = ^2 a Ji ®bV A ! 7 . (ii) 

We now apply this recipe to W m = V(t). Obviously, 
after a couple of products, W m potentially requires an ex- 
ponentially large auxiliary dimension to represent it per- 
fectly. It is here that we use a method similar to the 
DMRG truncation to reduce the size of this auxiliary di- 
mension. We begin by representing W l as an MPO per- 
fectly for as large an I as possible. Then we minimise 



the Hilbert-Schmidt norm difference \\W — V||hs = 
y/tT({W l - Yy{W l - Y)), where Y is a MPO with a 
smaller auxiliary dimension. This is a multiquadratic op- 
timisation problem and can be solved numerically in poly- 
nomial resources in n. (For a detailed description of this 
procedure see |7].) We then use the approximation Y to 
obtain an approximation YW to W l+1 and repeat this pro- 
cess for the desired number of iterations. 

Given an approximation V(f) = 
EjeQ„ A jo (t)A h (t) • • • A- ? "- 1 (i)cr- ? '°(8)cr- ? ' 1 (g>- • -^a^- 1 
to U(t) as an MPO with bounded auxiliary dimension 
D it is straightforward to obtain the trace efficiently: 
tr(a Q ) = 6 afi gives tr(V(t)) = A°A° • • • A° n _ v 

This procedure provides us with a discrete representation 
g k ~ ti(W k ), k = 0, 1, . . . , m — 1, for an approximation 
to the trace f{t k ) of U{t k ). To obtain an approximate rep- 
resentation for fi H we apply a discrete fourier transform to 

g k = e ' g k for some s chosen small enough to cutoff 
g k completely by /c max . A standard result of fourier anal- 
ysis shows that this discrete representation will be a good 
representation for the smeared eigenvalue density function 
as long as we sample more rapidly than the Nyquist fre- 
quency v. We can provide a bound for v by noting that 
the largest eigenvalue E max of H satisfies E max < \\H\\, 
which is 0(n). 

We have applied this numerical method to study the 
eigenvalue distribution function for a random antiferro- 
magnetic Heisenberg model on 20 spins, see Fig. |2] and 

Fig. 

Helpful discussions with Jens Eisert and Frank Ver- 
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